package 曲线多项式;
import java.util.Arrays;
import java.util.Vector;
import org.ejml.simple.SimpleMatrix;
//三次样条插值
public class cubic_spline_interpolation {
     // 样条函数有 4个参数  ai、bi、ci、di
     // Si(x) = ai+bi(x − xi)+ci(x − xi)(x − xi) +di(x − xi )(x − xi )(x − xi )
     private final Vector<Double> ai = new Vector<Double>(); // ai、bi、ci、di
     private final Vector<Double> bi = new Vector<Double>();
     private final Vector<Double> ci = new Vector<Double>();
     private final Vector<Double> di = new Vector<Double>();
     // 笛卡尔坐标的坐标点列表 x、y
     public Vector<Double> xi = new Vector<Double>();
     public Vector<Double> yi = new Vector<Double>();
     private int len;// 计算初始 坐标点的个数
     private Double[] toarray;
     public cubic_spline_interpolation(){}
     public cubic_spline_interpolation(double[] x, double[] y){
          this.len = x.length;
          this.toarray = new Double[this.len];
          Vector<Double> hi = new Vector<Double>();  // 计算相邻参考点间的 x 差值；
          for (int i = 0; i < this.len; i++) {
               this.xi.add(i, x[i]);                  //  求解 xi
               this.yi.add(i, y[i]);                  //  求解 yi
               this.ai.add(i, y[i]);                  //  求解 ai
               if (i < this.len - 1) {
                    hi.add(i, x[i + 1] - x[i]);
               }
          }
          SimpleMatrix A = calcA(hi);
          SimpleMatrix B = calcB(hi);
          SimpleMatrix temp = A.solve(B);
          for (int i = 0; i < this.len; i++) {
               this.ci.add(i,temp.get(i,0));      //  求解 ci
          }
          double tb;
          for (int i = 0; i <  this.len-1; i++) {     //  求解 di、bi
               this.di.add(i,(this.ci.get(i+1) - this.ci.get(i)) / (3.0*hi.get(i)) );
               tb = (this.ai.get(i+1) - this.ai.get(i)) / hi.get(i) - hi.get(i)*(this.ci.get(i+1) + 2*this.ci.get(i))/3.0;
               this.bi.add(i,tb);
          }
     }
     private SimpleMatrix calcA(Vector<Double> hi) {  // 计算A矩阵来求解三次样条的参数 ci
          SimpleMatrix A = new SimpleMatrix(this.len, this.len); // 初始化 len 行 len 列 的A矩阵
          A.set(0, 0, 1.0);
          for (int i = 0; i < this.len - 1; i++) {
               if (i != this.len - 2) {
                    A.set(i + 1, i + 1, 2 * (hi.get(i) + hi.get(i + 1)));
               }
               A.set(i + 1, i, hi.get(i));
               A.set(i, i + 1, hi.get(i));
          }
          A.set(0, 1, 0.0);
          A.set(this.len - 1, this.len - 2, 0.0);
          A.set(this.len - 1, this.len - 1, 1.0);
          return A;
     }
     private SimpleMatrix calcB(Vector<Double> hi){       // 计算B矩阵来求解三次样条的参数 ci
          SimpleMatrix B = new SimpleMatrix(this.len,1);
          for (int i = 0; i < this.len - 2; i++){
               B.set(i+1, 0,
                   3*(this.ai.get(i+2) - this.ai.get(i+1))/hi.get(i+1) -
                   3*(this.ai.get(i+1) - this.ai.get(i))/hi.get(i)
               );
          }
          return B;
     }
     private int findIndex(double x){  //找到 x 所属的 xi 的序列号
          this.xi.toArray(this.toarray);
          int index = Arrays.binarySearch(this.toarray, x);
          if (index < 0) {
               index = -(index + 2);// 如果未找到，计算插入点
               return index;
          }
          return index;
     }
     private boolean isBorder(double x){  //  判断越界
         return x < this.xi.get(0) || x > this.xi.lastElement();
     }
     public double calc(double x){  // 根据 x 拟合 y 值
          if ( isBorder(x) ){
               System.out.println("参数越界");
               return Double.POSITIVE_INFINITY;
          }
          int temp = findIndex(x);
          double dx = x - this.xi.get(temp);
          return this.ai.get(temp) + this.bi.get(temp)*dx + this.ci.get(temp)*Math.pow(dx,2) +  this.di.get(temp)*Math.pow(dx,3);  // Si(x) = ai+bi(x − xi)+ci(x − xi)(x − xi) +di(x − xi )(x − xi )(x − xi )
     }
     public double calcd(double x){  //  拟合 y 值的一阶导数
          if ( isBorder(x) ){
               System.out.println("参数越界");
               return Double.POSITIVE_INFINITY;
          }
          int temp = findIndex(x);
          double dx = x - this.xi.get(temp);
          return this.bi.get(temp) + 2*this.ci.get(temp)*dx + 3*this.di.get(temp)*Math.pow(dx,2);
     }

     public double calcdd(double x){  //  拟合 y 值的二阶导数
          if ( isBorder(x) ){
               System.out.println("参数越界");
               return Double.POSITIVE_INFINITY;
          }
          int temp = findIndex(x);
          double dx = x - this.xi.get(temp);
          return  2*this.ci.get(temp) + 6*this.di.get(temp)*dx;
     }

}